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O ■ ABSTRACT 



Hollow-core holey fibers are promising candidates for low- loss guidance of light in various applications, e.g., for 
Q the use in laser guide star adaptive optics systems in optical astronomy. We present an accurate and fast method 
^ _ for the computation of light modes in arbitrarily shaped waveguides. Maxwell's equations are discretized using 
Q . vectorial finite elements (FEM). We discuss how we utilize concepts like adaptive grid refinement, higher-order 
' ^ ' finite elements, and transparent boundary conditions for the computation of leaky modes in photonic crystal 
fibers. Further, we investigate the convergence behavior of our methods. 

We employ our FEM solver to design hollow-core photonic crystal fibers (HCPCF) whose cores are formed 
from 19 omitted cladding unit cells. We optimize the fiber geometry for minimal attenuation using multidimen- 
sional optimization taking into account radiation loss (leaky modes). 
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^ 1. INTRODUCTION 

Photonic crystal fibers (PCF) have a core surrounded by a periodic arrangement of holes and struts. This 
: structure prevents leakage of light to the exterior.^ In contrast to ordinary index- guiding optical fibers the core 
can have a smaller refractive index than the surrounding material or even be hollow. Although hard to fabricate 
1 there are many application areas where the guidance principles of PCFs offer great advantages. An important 
Q ■ example for the use of hollow-core photonic crystal fibers (HCPCF) is the field of high-power light transmission,^ 
, e.g. for pulsed lasers for adaptive optics systems in astronomy.^ Loss mechanisms of HCPCFs are coupling of 
J>-^ the fundamental mode to interface modes due to diffuse scattering at glass/air interfaces and leakage of light 
, to the exterior through the photonic crystal structure.^ Another problem is the usually narrow transmission 
bandwidth. 

^ When producing a fiber it is therefore very important to be able to optimize the fiber design to match desired 

properties and minimize undesired losses. Because of the complicated structure of PCFs this can only be done 
5J] numerically. The oldest approach for computation of modes in fibers is the plane- wave expansion (PWE) method. 
5^ , Because of the complicated structure and great number of refractive index jumps in the fiber cross section this 
method is inefficient. For example discontinuous behaviour of the propagation modes at the glass/air interfaces 
can only be described with a huge number of basis functions taken into account for computation. This leads to 
very large computation times and poor convergence behaviour of the method.^ Moreover, the modeling of an 
infinite exterior and therefore the simulation of radiation leakage is very difficult with the PWE method. Here 
we present the finite-element method (FEM) for the computation of leaky modes in photonic crystal fibers. 
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2. FORMULATION OF THE PROPAGATION MODE PROBLEM 



In this section wc will derive the mathematical formulation of the propagation mode problem. The geometry of 
a fiber is invariant in one spacial dimension (along the fiber), here z-direction. A propagating mode is a solution 
to the time harmonic Maxwell's equations, which exhibits a harmonic dependency in ^-direction: 



E = Epni(.x,j/)exp(2fc2z) 
H = llp„^{x,y)cxp{ik^z) . 



(1) 



Epm(a;,y) and Hpni(a;,y) are the electric and magnetic propagation modes and the parameter kz is called prop- 
agation constant. If the permittivity e and permeability fj, can be written as: 
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we can split the propagation mode into a transversal and longitudinal component: 

Hpm(a;,t/) 
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Inserting (1) with (2) and (3) into Maxwell's equations yields: 
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Now we define H^ = k^H^ and get: 
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Eq. (6) is a quadratic eigenvalue problem for the propagation constant and propagation mode Hpni(.T,y). 
We get a similar equation for Ep„i(a;,2/) exchanging e and /z. For our numerical analysis we will not look at the 
propagation constant but define the effective refractive index rieff which we will also refer to as eigenvalue: 
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(7) 



where Aq is the vacuum wavelength of light. 



3. SOLUTION OF THE PROPAGATION MODE PROBLEM WITH THE 

FINITE-ELEMENT METHOD 

In order to compute propagation modes we have to solve the eigenvalue problem (6) numerically. Note that 
this problem is formulated on R^. We have an eigenvalue problem on an unbounded domain. Usually, propa- 
gation modes are computed on an artificially bounded domain applying either periodic or so called Dirichlet 
boundary conditions, i.e. Epi„(x,2/) = for {x,y) G dSl (propagation mode vanishes on the boundary dCl of the 
computational domain). Using this simplification it is not possible to model leaking of light from the core of 
the fiber to the exterior (leaky modes). Here we want to take this effect into account. Since our computational 
domain still has to be of finite size, we apply so-called transparent boundary conditions to dCl. We realize these 
boundary conditions with the perfectly matched layer (PML) method.^ The propagation constant kz becomes 
complex and the corresponding mode is damped according to exp{—^{kz)z) while propagating along the fiber. 
This dampening is due to radiation leakage from the fiber to the exterior. 



(a) (b) 
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Figure 1. (a) Geometry of HCPCF used for mode computation; (b) boundary conditions for mode computation with a 
quarter of the fiber; (c) geometrical parameters describing HCPCF. Pitch A, hole edge radius r, strut thickness w, core 
surround thickness t; (d) detail from a triangulation of HCPCF. Due to the flexibility of triangulations all geometrical 
features of the HCPCF are resolved. 

In contrast to the PWE method, whose ansatz functions are spread over the whole computational domain 
(plane waves) the FEM method uses localized ansatz functions. To construct these ansatz functions, the compu- 
tational domain has to be discretized. This means the geometry is subdivided into N patches, e.g. triangles in 
two dimensions, tetrahedra in three dimensions. Fig. 1(d) shows such a triangulation of a photonic crystal fiber 
(Fig. 1(c)). On the patches usually polynomial ansatz functions are defined. Since we are solving Maxwell's 
equations, our solution is a vectorial function and therefore we use vectorial ansatz functions Vi(x^y). The 
numerical solution for the electric field E is a superposition of these localized ansatz functions 

N 

^{x,y) =^hvi{x,y) (8) 

1=1 

The FEM computation determines the unknown coefficients hi. The FEM method has several advantages: 

• Maxwell's equations are solved rigorously without approximations. 

• The flexibility of triangulations allows the computation of virtually arbitrary structures without simplifi- 
cations or approximations, as illustrated in Fig. 1(d). 

• Choosing appropriate ansatz functions Vi(x,y) for the solution of Maxwell's equations, physical properties 
of the electric field like discontinuities or singularities can be modeled very accurately and don't give rise 
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Figure 2. First, second and fourth fundamental core modes of HCPCF illustrated in Fig. 1(a) - Parameters: A = 1550 nm, 
r — 300 nm, w = 50 nm, t = 170 nm. 

to numerical problems. Such discontinuities often appear at glass/air interfaces of photonic crystal fibers, 
see Fig. 2(c). 

• Adaptive mesh-refinement strategies lead to very accurate results and small computational times which is 
a crucial point for optimization of fiber design. 

• The FEM approach converges with a fixed convergence rate towards the exact solution of Maxwell-type 
problems for decreasing mesh width (i.e. increasing number N of sub-patches) of the triangulation. There- 
fore, it is easy to check if numerical results can be trusted.^ 

In the following sections we apply the FEM method to the computation of leaky propagation modes in HCPCFs. 
The only simplification we make is to extend the fiber cladding to infinity and thereby neglect its finite size. 
This is justified if the cladding of the fiber is much larger than the microstructured core, and no light entering 
the cladding is reflected back into the core, which is usually the case. 

Throughout this paper we use the FEM package JCMsuite developed at the Zuse Institute Berlin for the 
numerical solution of Maxwell's equations. It has been successfully applied to a wide range of electromagnetic 
field computations including waveguide structures,® DUV phase masks, ^ and other nano-structured materials.^' ^ 
It provides higher-order edge elements, multi-grid methods, a-posteriori error control, adaptive mesh refinement, 
and adaptive transparent boundary conditions. Further numerical details about the computation of leaky modes 
with JCMsuite can be fomid in.^° 

4. NUMERICAL ASPECTS OF PROPAGATION MODE COMPUTATION 

Fig. 2 shows the first, second, and fourth core modes with lowest energy of the HCPCF illustrated in Fig. 1(a). 
All of these core modes appear several times in the spectrum because of the Cgv symmetry of the layout. To 
decrease the problem size, we only use one quarter of the fiber as computational domain. Of course we should 
obtain the same eigenvalues taking the full, half or quarter fiber as computational domain, as demonstrated in 
Table 1. The imaginary part of the eigenvalue is much smaller than the real part. The differences in 3(nofr) are 
therefore due to the numerical uncertainty at the chosen refinement level. 3(7ioff) will of course converge with 
increasing number of ansatz functions, see Fig. 3. Figure 1(b) shows which boundary conditions we apply to 
the artificial inner boundaries for the case of a quarter fiber. These boundary conditions follow from symmetry 
of the core modes. Note that the imaginary and real parts of the eigenvalues differ by 11 orders of magnitude. 
Computation of leaky modes is therefore a multi-scale problem which makes it numerically very difficult. When 
analyzing radiation losses of a fiber, the imaginary part of the effective refractive index riog is the quantity of 
interest. In the following we analyze how accurately we can compute Ucs, i.e. we look at the convergence of the 
FEM computation. Therefore we compute UeS with an increasing number of unknowns. We can either refine 





unknowns 


1st eigenvalue 


2nd eigenvalue 


3rd eigenvalue 


full fiber 


861289 


0.998265726+ 9.508 -lO-^^i 


0.99212759+ 1.513- 10-^"i 


0.9908968 + 2.089 •lO-^'Ji 


half fiber 


438297 


0.998265719+ 9.357- 10-^^i 


0.99212742+ 1.489- lO-i'^i 


0.9908958 + 2.086 -lO-^^i 


quarter fiber 


218504 


0.998265724+ 9.372 • IQ-'^'H 


0.99212652+ 1.504- 10"^"i 


0.9909169 + 2.302 -10"^"i 



Table 1. First, second and third eigenvalue computed with full, half and quarter fiber as computational domain. 
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Figure 3. Relative error of first eigenmode in dependence on number of unknowns of FEM computation for difi^erent FEM 
degrees p. (a) and (b): adaptive refinement; (c) and (d) comparison of adaptive and uniform refinement. Parameters: 
A = 1550 nm, r = 300 nm, w = 50 nm, t = 170 nm. 



the triangulation or increase the polynomial degree p of the ansatz functions on each patch to obtain a more 
accurate numerical result. The software package JCMsuite offers finite elements with a maximum degree of 9. 
Furthermore, the refinement of the grid can be performed uniformly or adaptively. In a uniform refinement step 
each triangle is subdivided into 4 smaller ones. In an adaptive refinement step, an a-posteriori error estimator 
automatically chooses a certain number of triangles where the error of the FEM solution is large. Only these 
triangles are refined. This method leads to accurate results with smaller number of unknowns which also implies 
shorter computational times and smaller memory requirements. Fig. 3 shows the relative error of the real 
and imaginary parts of rics in dependence on the number of unknowns of the FEM computation for different 
FEM degrees p . The relative error decreases with an increasing number of unknowns. Let us first look at the 



convergence with adaptive refinement 3(a), (b). For the real part we achieve the fastest convergence for a finite 
element degree of 2. Note that already for a finite element degree of 1 and only 20, 000 unknowns we have a 
relative error smaller than 10~^. For 10^ unknowns we have a relative error smaller than 10~® for FEM degrees 
greater than 1. The imaginary part of the effective refractive index is much harder to compute accurately. For 
the same number of unknowns its relative error is much larger than the relative error of the real part. For a 
FEM degree of 1 the imaginary part does not converge at all up to 3 • 10^ unknowns. However the higher the 
FEM degree the faster the imaginary part converges. For the computation of leaky modes in HCPCFs it is 
therefore important to employ higher-order finte elements. Fig. 3(c), (d) shows a comparison of the convergence 
of the real and imaginary parts between adaptive and uniform refinement. While adaptive refinement leads to 
a faster convergence of the real part for all FEM degrees the imaginary part converges almost equally fast for 
an FEM degree of 4 and even slower for an FEM degree of 2. This happens because the imaginary part of 
the eigenvalue is much smaller than the real part. The error estimator which chooses triangles in the case of 
adaptive refinement primarily reduces the error of the real part first, since its contribution to the relative error 
of the complex eigenvalue is much larger than the imaginary part. Computing leaky modes a uniform refinement 
strategy is therefore more suitable if one is interested in the imaginary part and radiation losses. 

For further numerical analysis of radiation losses of the HCPCF we choose a finite-element degree of 6 and 
no refinement step. With these settings the numerical problem has approximately 500, 000 unknowns for a 
triangulation according to Fig. 1(d). The relative error of the real and imaginary parts is « 10~^ and ~ 10~^ 
respectively, according to the convergence curves. Such a computation requires 2GB RAM and about 10 minutes 
computation time on a standard desktop PC. It is fast enough for multidimensional optimization of the fiber 
design. 

5. OPTIMIZATION OF HCPCF DESIGN 

In this section we want to optimize the fiber design shown in figure 1(a) in order to reduce radiation losses, 
i.e. we want to minimize Sj(neff). The basic fiber layout is a 19-cell core with rings of hexagonal cladding cells. 
Since these cladding rings prevent leakage of radiation to the exterior we expect that an increasing number of 
cladding rings reduces radiation leakage and therefore Q(nefr). This is confirmed by our numerical simulations 
shown in Fig. 4. The radiation leakage decreases exponentially with the number of cladding rings and thereby 
the thickness of the photonic crystal structure. This behavior agrees with the exponential dampening of light 
propagating through a photonic crystal structure with frequency in the photonic band gap. 

For our further analysis we fix the number of cladding rings to 6. The free geometrical parameters arc the 
pitch A, hole edge radius r, strut thickness w, and core surround thickness t, see Fig. 1(c). Fig. 4 shows 
the imaginary part of the effective refractive index in dependence on these parameters. For each scan all but 
one parameter were fixed. For the strut thickness w and the hole edge radius r we find well-defined optimal 
values which minimize SJ(neff)- For pitch A and core surround thickness t a large number of local minima and 
maxima can be seen. Now we want to optimize the fiber design using multidimensional optimization with the 
Nelder-Mead simplex method.^ To reduce the number of optimization parameters we fix the hole edge radius to 
the determined minimum at r = 354 nm since its variation has the smallest effect on $5(neff). For optimization 
we have to choose starting values for A, t and w. Since the simplex method searches for local minima we have 
to decide in which local minimum of A and t we want to search. We choose t = 152 nm since here 3(noff) has 
a global minimum and A = 1550 nm since the bandwidth of this minimum is much larger than for the global 
minimum at A = 1700 nm. 

Optimization yields a minimum value of 3(noff) = 5 • 10~^^^ for the imaginary part of the effective refractive 
index. The corresponding geometrical parameters are A = 1597 nm, w = 38 nm , t = 151 nm. 

6. CONCLUSION 

We have demonstrated that the finite-element method is very well suited for the analysis of light propagation in 

hollow-core photonic crystal fibers. Here we considered a 19-ccll HCPCF with 6 cladding rings. A convergence 
analysis allowed us to quantify the relative error of real and imaginary parts of the propagation constant of leaky 
eigenmodes. The real part could be computed in 20 s on a standard desktop PC down to a relative error of 
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Figure 4. Imaginary part of effective refractive index 9(neff ) in dependence on: (a) number of cladding rings, (b) pitch 
A, (c) core surround thickness t, (d) strut thickness w, (e) hole edge radius r. Parameters: A = 1550 nm, r = 300 nm, 
w = 50 nm, t = 170 nm, 6 cladding rings, wavelength A = 589 nm. 

10"''. The imaginary part, being 10 orders of magnitude smaller, could be determined in about lOmin with a 
relative error smaller than 10""^. Thanks to the short computational time, we were able to optimize the fiber 
design regarding pitch, strut thickness, core surround thickness, and hole edge radius of the cladding cells and 
therewith minimize radiation losses. 
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